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I.  INTRODUCTION 


Deterministic  simulation  of  earthquake  ground  motion  has 
played  an  Increasingly  Important  role  in  seismology  and  earthquake 
engineering  In  recent  years.  For  example,  ground  motion  simulation 
has  been  used  recently  as  a  tool  for  developing  engineering  design 
motion  criteria  (Wiggins,  et  aU ,  1978;  Apsel,  1979).  Such  simula¬ 
tions  require  theoretical  models  for  both  the  source  process  and  the 
propagation  and  dissipation  of  seismic  energy.  While  ground  motion 
simulations  have  been  undertaken  using  rather  rigorous  theoretical 
methods  to  model  anelastlc  wave  propagation  from  source  to  site, 
including  the  effects  of  depth-dependent  geologic  structure,  the 
earthquake  source  Itself  Is  usually  specified  on  largely  Intuitive 
grounds.  The  displacement-discontinuity  time  history  (slip  func¬ 
tion),  by  means  of  which  the  earthquake  source  Is  represented.  Is 
generally  prescribed  without  rigorous  consideration  of  fault 
dynamics. 

Following  Haskell  (1964;  1969)  and  Savage  (1966),  most  studies 
have  represented  the  earthquake  by  a  slip  function  which  is 
spatially  uniform  over  the  fault  surface  and  has  some  simple  time 
dependence,  usually  that  of  a  finite-duration  ramp.  Simple  kine¬ 
matic  source  models  of  this  type  have  proven  useful  for  representing 
low-frequency  (less  than  about  1  Hz)  characteristics  of  earthquake 
ground  motion.  For  example,  Bouchon  (1979)  used  a  uniform  disloca¬ 
tion  of  ramp  form  in  combi  nation  with  a  layered  earth  structure  to 
model  the  Station  2  velocity  and  displacement  recordings  of  the  1966 
Parkfleld  earthquake.  However,  high  frequency  (1  to  20  Hz)  ground 
motion  Is  highly  sensitive  to  the  specification  of  the  source 
process.  The  analysis  of  Madriaga  (1978),  in  particular, 
underscores  the  inapprorlateness  of  the  uniform-slip  kinematic 
models  for  synthesizing  ground  motion  with  wavelengths  much  shorter 
than  the  fault  widtn.  This  high  frequency  radiation  is  important  in 
earthquake  engineering,  as  well  as  in  nuclear  monitoring  studies. 


1 


SYSTEMS.  SCIENCE  AND  SOETWAEE 


r 


Closed-form  theoretical  solutions  describing  the  slip  function 
are  available  only  for  the  most  Idealized  dynamic  shear-crack  models 
of  earthquakes.  Restricting  consideration  to  three-dimensional 
analyses,  perhaps  the  most  useful  of  such  analytical  results  are  the 
solutions  of  Kostrov  (1964)  and  Burrldge  and  Willis  (1969).  The 
former  gives  the  slip  history  on  a  circular  shear  crack  which 
initiates  at  a  point  In  a  prestressed  whole  space  and  grows  at  a 
fixed  rupture  velocity  without  stopping;  the  latter  extends 
Kostrov's  result  to  the  case  of  an  elliptical  shear  crack.  These 
self-similar  solutions  are  characterized  by  square-root  singular¬ 
ities  at  the  crack  edge  In  both  shear  stress  and  slip  velocity,  and 
the  intensity  of  these  singularities  grows  as  the  crack  dimension 
Increases.  While  these  analytic  solutions  are  very  useful  for 
interpreting  tne  results  of  more  complex  numerical  studies,  they 
cannot  account  for  effects  associated  with  the  stopping  of  rupture 
growth  and  the  ensuing  arrest  of  slip. 

Madariaga  (1977a)  used  two-dimensional  analytical  results, 
notably  those  of  Kostrov  (1966,  1975)  and  Possum  and  Freund  (1975) 
to  characterize  the  slip-velocity  singularities  associated  with  the 
starting  and  stopping  of  ruptures.  He  then  applied  a  representation 
theorem,  together  with  Keller's  (1962)  geometric  theory  of  dif¬ 
fraction  to  construct  a  high-frequency  approximation  for  the  radi¬ 
ation  from  shear  cracks  In  three  dimensions.  This  analysis  yields 
expressions  for  the  radiation  from  a  discrete  jump  in  rupture  velo¬ 
city.  The  solution  Involves  a  stress  intensity  factor  which  depends 
on  the  three-dimensional  geometry  of  the  fault  as  well  as  its  stress 
and  rupture  history,  and  a  second  factor  which  depends  only  on  the 
instantaneous  jump  in  rupture  velocity  along  the  crack  edge.  The 
solution  In  this  form  provides  considerable  Insight  Into  the  process 
of  high-frequency  generation,  although  to  fully  characterize  the 
stress  Intensity  factor  and  rupture  velocity,  a  complete  solution  to 
the  three-dimensional  dynamical  problem  would  still  be  required. 


In  general,  numerical  methods  are  necessary  to  solve  the 
three-dimensional  dynamic  problem  of  a  fault  which  stops.  Several 
studies  have  addressed  this  problem,  with  the  approximation  that 
rupture  velocity  Is  specified  a  priori.  This  fixed-rupture-velocity 
fault  model  has  been  studied  for  faulting  confined  to  a  circular 
region  (Madariaga,  1976;  Oas,  1980),  a  semi-circular  region 
(Archuleta  and  Frazier,  1978),  and  rectangular  regions  (Madariaga, 
1977b,  1979;  Day,  1979;  Archuleta  and  Day,  1980).  Oas  (1981)  and 
Day  (1979)  have  studied  rectangular  faults  In  which  rupture  velocity 
Is  derived  from  a  fracture  criterion  (spontaneous  rupture).  These 
numerical  solutions  demonstrate  that  edge  effects  associated  with 
the  narrow  dimension  of  the  fault  substantially  Influence  the  slip 
function,  controlling  static  slip  and  slip  rise  time. 

In  this  study,  our  objective  is  to  provide  an  improved  under¬ 
standing  of  three-dimensional  geometrical  effects  governing  slip 
functions.  The  study  employs  three-dimensional  finite-difference 
solutions  to  the  dynamic  fault  problem  In  a  whole  space.  On  the 
basis  of  the  numerical  solutions,  closed-form  approximations  are 
derived  for  the  static  slip  and  rise  time  predicted  for  rectangular 
faults.  These  quantities  are  measures  of  the  low-  to  Intermediate- 
frequency  content  of  the  slip  function.  A  particular  emphasis  of 
this  study,  however,  will  be  to  quantify  the  high-frequency  behavior 
of  the  slip  function,  which  has  been  largely  unresol vable  from  pre¬ 
vious  numerical  studies  of  earthquake  dynamics.  It  Is  the  high- 
frequency  character  (greater  than  1  Hz  or  so)  of  the  slip  function 
which  Is  of  primary  Importance  for  synthesizing  near-field  ground 
motion  In  the  period  range  of  engineering  interest.  The  high-fre¬ 
quencies  are  also  of  Importance  for  synthesizing  earthquake  ground 
motion  at  regional  distances,  where  substantial  seismic  energy  is 
recorded  in  the  1-5  Hz  range.  A  good  source  model  for  high 
frequencies  Is  also  Important  for  inferring  the  attentuative 
properties  of  teleseismic  propagation  paths  (i.e.,  t*). 
Furthermore,  some  of  the  methods  proposed  for  discriminating 
earthquakes  from  underground  explosions  use  spectral  characteristics 
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of  teleslsmlc  P-waves  in  the  1-5  Hz  range  (e.g.,  Bache,  et  al., 
1979).  Once  the  high-frequency  behavior  of  the  slip  function  has 
been  quantified  from  the  numerical  solutions,  the  numerical  results 
are  combined  with  asymptotic  results  for  dynamic  cracks,  in  order  to 
characterize  the  stress  Intensity  and  energy  release  rate  at  the 
advancing  fault  edges. 

So  as  to  focus  our  analysis  on  the  geometrical  effects 
associated  with  the  three-dimensionality  of  the  problem,  we  employ 
the  approximation  that  rupture  velocity  be  a  specified  constant. 
The  finite  difference  method  used  here  is,  of  course,  also  well 
suited  to  the  more  complex  problem  of  modeling  spontaneous  rupture 
propagation. 


The  finite  difference  computations  on  which  this  analysis  is 
based  were  performed  on  the  ILLIAC  IV  computer  at  NASA/Ames  Research 
Center.  Programming  support  was  provided  under  the  direction  of 
Susan  Blester  and  Stewart  Hopkins. 
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II.  THE  FAULT  MODEL 


2.1  INTRODUCTION 

In  the  earthquake  simulations  reported  here,  we  treat  faulting 
as  a  propagating  shear  stress  relaxation  which  occurs  as  a  conse¬ 
quence  of  shear  failure  on  a  planar  surface.  The  mathematical 
formulation  follows  closely  that  of  Equations  2.1  to  2.13  of  Kostrov 
(1970).  Archuleta  and  Frazier  (1978)  also  present  a  detailed 
exposition  of  a  mathematical  model  of  a  propagating  shear  stress 
relaxation. 

We  will  specify  the  Initial  state  of  stress  in  the  medium,  its 
constitutive  properties,  the  rupture  velocity,  the  limiting  edge  of 
the  rupture  surface,  and  the  friction  law  to  be  satisfied  on  the 
rupture  surface  after  failure.  Although  the  velocity  of  the  rupture 

front  Is  prescribed,  the  time  of  arres£  of  slip  at  a  given  point  is 
not.  Instead,  the  cessation  of  slip  is  a  consequence  of  the 
nonlinear  friction  law,  and  is  determined  as  a  part  of  the  dynamic 
solution. 

2.2  INITIAL  CONDITIONS  AND  CONSTITUTIVE  PROPERTIES 

For  time  t  less  than  zero,  we  assume  that  an  equilibrium  state 
of  stress  exists  with  velocity  everywhere  zero.  The  equilibrium 
configuration  is  such  that  the  prospective  fault  plane  experiences  a 
uniform  shear  traction  and  compresslonal  normal  traction  <ru. 

The  fault  plane  Is  permitted  to  fail  in  shear,  but  the  medium 
will  otherwise  be  assumed  to  be  linearly  elastic.  Since  average 
stress  changes  associated  with  faulting  are  modest,  on  the  order  of 
a  few  hundred  bars,  linear  elasticity  is  a  reasonable  model  of 
material  behavior  away  from  the  immediate  zone  of  faulting. 

2.3  GROWTH  OF  THE  RUPTURE 

The  rupture  surface  is  assumed  to  occupy  a  prescribed  plane 
with  unit  normal  vector  n.  We  specify  the  growth  of  the  fault 
surface  as  a  function  of  time,  rather  than  determining  its  evolution 
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from  the  dynamic  solution  via  some  failure  model.  The  rupture 
nucleates  at  a  point  and  expands  symmetrically  at  a  constant, 
prescribed  rupture  velocity  u^.  until  It  reaches  a  prescribed 
rectangular  boundary  (Figure  1).  E(t)  denotes  the  portion  of  the 
plane  which  has  ruptured  by  time  t;  w  and  i  denote  the  width  and 
length  of  a  rectangular  fault,  and  x  and  y  are  Cartesian  coordinates 
on  the  fault  plane.  Then  E(t)  consists  of  all  points  x,  y  such  that: 


lx! 


t 


and  |y|  s<  J  . 

Two  somewhat  artificial  features  of  this  rupture  model  require 
mention.  They  are:  1)  the  instantaneous  acceleration  of  the  rupture 
to  its  terminal  velocity,  and  (2)  the  Instantaneous  deceleration  of 
rupture  velocity  to  zero  along  a  prescribed  boundary.  The  former 
assumption  may  be  a  fairly  good  approximation.  There  is  both 
experimental  (Mu,  et  al_. ,  1972:  Archuleta  and  Brune,  1975)  and 
theoretical  (Cherry,  et  al_.,  1976;  Das  and  Aki,  1977)  evidence  that 
rupture  velocity  can  accelerate  very  rapidly  to  its  terminal  value. 

The  approximation  of  abrupt  stopping,  on  the  other  hand,  is 
difficult  to  support  experimentally,  since  ruptures  normally 
propagate  completely  through  laboratory  samples.  While  the 
approximation  of  abrupt  stopping  may  be  quite  artificial,  it  Is  not 
precluded  by  theories  of  dynamic  crack  propagation.  For  example, 
Husselni  et  al_. ,  (1975)  have  shown  that  a  rupture  can  stop 

instantaneously  when  it  encounters  jumps  In  fracture  energy  on  the 
fault  plane.  This  reflects  the  fact  that  a  crack  edge,  at  least  in 
the  linearly-elastlc  continuum  theory,  lacks  "inertia".  That  is, 
the  stresses  imnediately  ahead  of  the  crack  edge  depend  on  rupture 
velocity,  but  not  on  the  time  derivatives  of  rupture  velocity 
(Eshelby,  1969). 
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Figure  1.  Rupture  geometry  and  coordinate  system  for  the 
numerical  simulations.  The  shear  prestress  Is  In  the  x 
direction  on  the  plane  z  *  0.  Rupture  Initiates  at  the 
origin  and  expands  symmetrically  at  fixed  rupture 
velocity* 
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The  Importance  of  the  mode  of  stopping  lies  In  Its 

consequences  for  high-frequency  radiation.  Madariaga  (1977a)  has 

shown  that  the  strongest  radiation  of  high-frequencies  Is  associated 

with  abrupt  changes  In  rupture  velocity  such  as  the  sudden  stopping 

at  the  fault  edges  In  our  Model.  A  rupture  velocity  jump  generates 

.2 

f  high-frequency  behavior  of  the  far-fleld  displacement 

spectrum.  In  contrast  to  the  starting  phase,  which  (assuming 

-3 

nucleatlon  at  a  point)  generates  at  most  an  f  spectral  asymptote. 

2.4  BOUNDARY  CONDITION  OF  THE  FAULT 

On  E(t),  we  permit  a  tangential  displacement  discontinuity 
(slip)  s(x,t),  and  require  continuity  of  the  traction  vector  and  of 
the  normal  component  of  displacement.  The  shear  traction  on  E  obeys 
a  simple  Coulomb  friction  law.  The  physical  requirements  of  this 
friction  law  are:  1)  the  magnitude  of  the  shear  traction  on  E  is 
bounded  by  a  prescribed  sliding  friction  value  which  depends  only  on 
the  normal  traction,  and,  ii)  the  shear  traction  is  equal  In 
magnitude  to  the  sliding  friction  value  and  opposite  in  direction  to 
the  slip  velocity  vector  whenever  the  latter  is  non- zero. 

The  vector  i  denotes  the  shear  traction  exerted  on  the  positive 
side  of  E  by  the  negative  side  (where  the  direction  of  n  is  from 
the  negative  side  of  E  toward  the  positive).  We  define  x^.  to  be  a 
sliding  frictional  traction  whose  direction  opposes  the 
Instantaneous  slip  velocity  and  whose  amplitude  is  proportional  to 
the  normal  traction  on  E  : 


where  x^  is  -u^,  the  product  of  the  normal  traction  and  the 
coefficient  of  dynamic  friction.  The  amplitude  of  the  sliding 
friction,  x^,  Is  presumed  to  be  positive,  constant,  and  less  than 
the  absolute  value  of  the  shear  prestress  x  ,  so  that  a  stress 
drop  occurs  at  the  rupture  front.  We  then  define  £c  to  be  the 
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shear  traction,  at  a  point  on  I,  which  would  be  sufficient  to 
enforce  continuity  of  velocity.  That  Is,  t.  Is  the  Instantaneous 
shear  traction  which  would  accompany  healing,  a  quantity  which  can 
be  readily  determined  at  any  time  step  from  the  numerical  solution. 
Then  the  following  boundary  condition  on  Z  Is  equivalent  to  the 
friction  law  described  In  the  last  paragraph: 

T.Uf  "  Hcl>Tf  '  (1) 

-  lie  ff  It,.]  «  Tf  . 

Equation  (1)  ensures  that  the  slip  velocity  £  Is  non-zero  only  If 
the  magnitude  of  t  would  otherwise  exceed  that  of  sliding  friction, 
xf. 


Note  from  Equation  (1)  that  we  do  not  assume  that  the 
frictional  stress  immediately  increases  above  the  dynamic  friction 
value  rf  when  a  slip  velocity  zero  occurs.  Rather  than  assuming 
instantaneous  recovery  of  strength,  we  let  rf  continue  to  bound 
the  shear  stress  on  the  fault,  even  after  the  slip  velocity  goes  to 
zero.  This  Is  In  accord  with  results  from  laboratory  studies  of 
time-dependent  rock  friction.  For  example,  Dieterich  (1972,  19781 
found  frictional  strength  of  rock  surfaces  to  be  proportional  to  the 
logarithm  of  their  time  of  contact,  with  essentially  no  recovery  of 
strength  occurlng  during  the  first  one  second  of  contact.  Thus,  any 
increase  of  frictional  strength  due  to  stationary  contact  should  be 
negligible  on  the  time  scale  of  dynamic  rupture. 

2.5  STRESS-DROP  SCALING 

The  dynamic  stress-drop,  ax,  is  defined  here  as  the  difference 
between  the  absolute  values  of  shear  prestress  and  sliding 
frictional  stress, 

ax  *  x0  -  Tf  . 

This  quantity  Is  the  stress  available  to  accelerate  the  fault  slip, 
and  has  also  been  termed  "effective  stress"  in  the  seismologies! 
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literature  (e.g.,  Brune,  1970).  If  the  slip  direction  were 
constrained  to  be  parallel  to  the  direction  of  the  prestress,  then 
the  dynamic  solution  would  scale  directly  with  the  assumed  value  of 
at,  and  would  be  independent  of  Tf.  This  constraint  of  the  slip 
direction  is  equivalent  to  assuming  t^»at.  In  the  current  study, 
we  have  made  this  simplifying  assumption  so  that  the  numerical 
results  can  be  scaled  rigorously  with  at.  The  assumption  t^»at 
is  consistent  with  laboratory  studies  of  stick-slip,  which  report 
fractional  stress  drops,  at/tq,  of  a  few  percent  to  a  few  tens  of 
percent  (for  example,  Byerlee,  1967;  Scholz,  et  aK,  1972; 
Dleterich,  et  al_. ,  1978). 

Actually,  this  simplification  does  not  have  a  significant 
effect  on  the  solution  (apart  from  inhibiting  slip  reversal  as 
discussed  below).  It  is  known,  for  example,  that  the  self-similar, 
expanding  elliptical  crack  exhibits  slip  which  is  everywhere 
parallel  to  the  prestress  direction,  even  for  the  case  of  zero 
friction  (Burridge  and  Willis,  1969).  For  finite  faults  with  low 
friction,  stopping  of  rupture  can  introduce  a  component  of  slip 
perpendicular  to  the  prestress  direction.  However,  Madariaga  (1976) 
has  shown,  for  the  case  of  a  finite  circular  crack,  that  this 
component  is  quite  small.  So,  scaling  with  at,  while  rigorous  only 
for  Tf»AT,  should  be  a  good  approximation  even  for  relatively  low 
values  of  Tf. 

2.6  HEALING 

When  the  slip  velocity  at  a  point  goes  to  zero.  Equation  (1) 
provides  the  criterion  for  whether  to  permit  further  slip. 
Recommencement  of  slip  is  prohibited  by  Equation  (1)  if  such  slip 
(which  must  be  accompanied  by  a  shear  traction  of  magnitude  t^) 
would  Increase  the  magnitude  of  the  shear  traction  x_,  rather  than 
decrease  it  since  that  would  violate  our  physical  assumption  that 
the  shear  traction  on  Z  opposes  the  slip  velocity. 
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Although  the  fraction  law  embodied  In  Equation  (1)  does  not 
automatically  preclude  a  reversal  of  sliding  direction,  a  reversal 
could  only  occur  for  very  small  values  of  sliding  friction  In  this 
model.  We  define  the  overshoot  stress  as  the  amount  by  which  the 
shear  traction,  at  a  point  on  Z,  drops  below  after  sliding 

ceases.  In  order  for  sliding  to  reverse,  the  shear  traction  would 
have  to  drop  from  the  sliding  friction  value  to  zero,  then  reverse 
sign  and  Increase  In  magnitude  bade  to  so  no  sliding  reversal 
will  occur  unless  Is  less  than  half  the  overshoot  stress.  The 
numerical  results  given  In  the  next  section  show  that  the  overshoot 
does  not  exceed  about  26  percent  of  the  dynamic  stress  drop; 
therefore,  a  sliding  reversal  would  only  occur  If  the  fractional 
stress  drop  at/tQ  were  at  least  0.88.  Since  we  have  assumed 
Tf»at  (i.e.,  very  small  fractional  stress  drop),  sliding  reversal 
does  not  occur  In  our  simulations. 
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III.  NUMERICAL  SOLUTIONS  FOR  THE  FAULT  SLIP 


3.1  METHOD  OF  SOLUTION 

The  mathematical  model  of  faulting  outlined  in  the  last 
section  poses  a  three-dimensional,  nonlinear,  mixed  boundary  value 
problem.  To  determine  the  fault  slip,  this  problem  Is  solved 
numerically  using  a  three-dimensional  finite  difference  method 
developed  by  Cherry  (1977),  In  which  explicit  time  stepping  Is  used 
to  integrate  the  dynamic  solution  in  time.  Equation  1,  governing 
the  fault  plane,  was  discretized  in  accordance  with  the  method 
developed  by  Oay  (1977,  Appendix  IV). 

3.2  DESCRIPTION  OF  THE  CALCULATIONS 

We  present  slip  functions  obtained  from  numerical  simulations 
of  four  different  fault  geometries.  Each  calculation  simulates  a 
rectangular  fault  surface  in  a  uniform  whole  space.  In  each  case, 
rupture  initiates  at  the  center  of  the  rectangle,  as  in  Figure  1, 
and  the  prestress  direction  is  aligned  with  the  long  dimension  of 
the  fault. 

The  numerical  results  have  been  scaled  to  represent  the 
following  set  of  physical  parameters: 

P  wave  speed,  a  *  6.0  km/sec 
S  wave  speed,  8  ■  3.46  km/sec 
shear  modulus,  u  *  3.24  x  10"  dynes/cm 
rupture  velocity,  uR  ■  3.12  km/sec 
dynamic  stress  drop,  at  ■  100  bars 

Three  of  the  calculations  represent  a  fault  length  £  of  8  km  and 
fault  widths  w  of  1.5,  4,  and  8  km,  respectively.  The  fourth 
calculation  was  for  a  fault  length  of  16  km  and  a  fault  width  of  4 
km.  The  solutions  may  be  rescaled  to  represent  a  different  set  of 
material  and  fault  parameters,  provided  Poisson's  ratio,  the  fault 
aspect  ratio,  and  the  ratio  uR/a  are  unchanged.  The  appropriate 
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scaling  relationships  have  been  summarized  by  Madariaga  (1976),  and 
we  repeat  then  here  for  convenience.  Assune  that  the  above  values 
of  dynamic  stress  drop,  fault  width,  shear  speed,  and  shear  modulus 
are  multiplied  by  factors  at',  w',  a',  and  u',  respectively.  Then, 
time,  length,  stress  and  dlsplacment  In  the  numerical  solutions  are 


to  be  multiplied  by 

scale  factors  t’ ,  x‘ ,  o'  and  u',  respectively: 

(time) 

t' 

.  JflL' 

6' 

(Length) 

x‘ 

-  w‘ 

(Stress) 

i 

0 

-  AT* 

(Displacement) 

U* 

II 

t 

• 

In  rescal 1 ng 

the  numerical  solutions,  however,  consideration 

must  be  given  to  a  hidden  length  scale  In  the  problem,  the  finite 
difference  cell  dimension.  The  main  limitation  of  the  finite 
difference  method  Is  that  the  discretization  causes  substantial 
numerical  dispersion  of  wavelengths  shorter  than  about  6  cell 
dimensions.  Consequently,  accuracy  is  degraded  for  high-frequency 
components  of  the  solution.  For  these  four  finite  difference 
calculations,  the  mesh  refinement  was  sufficient  to  retain  accuracy 
for  frequency  components  up  to  about  5  Hz.  Frequencies  higher  than 
5  Hz  have  therefore  been  removed  from  the  solutions  using  a 
'combination  of  artificial  viscosity  and  post-processing  with  a 
low-pass  digital  filter.  This  non-physical  cutoff  frequency  scales 
as  a/w. 

3.3  SLIP  TIME-HISTORIES 

Figure  2  shows  the  calculated  slip  histories  for  selected 
points  along  the  fault  length  for  the  square  (8  km  by  8  km)  fault. 
The  final  offset  Is  greatest  near  the  center  of  the  fault, 
decreasing  as  the  observation  point  approaches  an  edge.  The  rise 
time  of  the  slip  function  is  also  greatest  at  the  center,  where  It 
equals  approximately  the  fault  length  divided  by  the  shear  speed. 
The  rise  time  decreases  as  the  edges  of  the  fault  are  approached. 
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Computed  si fp  time-histories  along  the  center-line  of 
the  square  fault  plane  at  several  distances  from  the 
hypocenter.  Slip  Is  scaled  to  represent  a  dynamic 
stress  drop  of  100  bars  and  a  fault  length  of  8  km. 
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Arrest  of  slip  occurs  first  at  the  outer  edge  and  propogates 
Inward.  The  center  of  the  fault  is  the  last  point  to  heal. 

Overall,  the  healing  behavior  is  very  similar  to  that  of  the 

circular  fault,  as  obtained  by  Madariaga  (1976),  in  spite  of  the 

reduced  symmetry  present  In  the  square  fault  problem  (the  circular 
fault  stops  simultaneously  everywhere  on  its  perimeter  whereas  the 
square  fault  does  not). 

We  can  compare  the  residual  slip  for  the  square  fault  to 

Meuber's  (1937)  static  solution  for  a  circular  shear  crack  in  a 
Poisson  solid,  for  which  the  slip  s  is 

S(r)  -  a  (l-r2/a2)1/2  ,  (2) 

where  a  Is  the  crack  radius  and  r  Is  the  distance  to  the  center  of 
■;Me  crack.  For  the  circular  crack,  the  average  static  slip,  from 
Equation  2,  is  2/3  times  the  static  slip  at  the  center.  For  the 
re  fault  calculation,  the  average  static  slip  is  found  to  be 
0.65  times  the  static  slip  at  the  center,  which  is  nearly  identical 
to  the  circular  fault  relationship. 

As  a  reasonable  approximation,  we  might  apply  Equation  2  to 
the  square  fault  with  the  reinterpretation  that  a  is  V A /»’  ,  where 

A  is  the  fault  area.  In  that  case,  the  static  offset  at  the  center 
of  the  square  fault  exceeds  the  prediction  of  Equation  2  by  a  factor 
of  1.26.  This  "overshoot"  of  the  dynamic  solution  relative  to  the 
static  solution  has  been  discussed  for  circular  faults  by  Madariaga 
(1976),  Archuleta  (1976),  and  Das  (1980).  The  value  cf  1.26 
obtained  here  for  the  square  fault  Is  In  good  agreement  with  their 
numerical  results  (which  range  from  1.20  to  1.27). 

This  value  of  overshoot  implies  that  no  reversal  in  slip 
direction  will  occur  If  the  dynamic  frictional  traction exceeds 
about  0.13  at,  or  half  the  overshoot  stress.  Thus,  as  remarked  in 
the  previous  section,  a  slip  reversal  requires  a  fractional  stress 
drop,  at/tq,  greater  than  0.88,  which  Is  far  greater  than  those 
observed  in  laboratory  stick-slip  experiments.  We  emphasize  again 
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that  this  conclusion  requires  only  our  assumption  that  sliding 
friction  Is  a  dissipative  process  which  opposes  the  slip;*  It  does 
not  require  an  assumption  of  rapid  strength  recovery  on  the  fault 
following  a  slip  velocity  zero. 

Figure  3  shows  slip  histories  calculated  along  the  fault 
length  for  the  case  w  »  4  ,  £  ■  8.  The  time  histories  In  this  case 
are  similar  to  those  In  Figure  2,  but  the  rise  time  and  static  slip 
are  reduced,  especially  near  the  center  of  the  fault.  Rise  time  and 
static  slip  at  a  point  are  apparently  controlled  by  proximity  of  the 
nearest  edge  of  the  fault. 

Figure  4  shows  slip  histories  for  the  case  w  ■  1.5,  l  *  8. 
Again,  It  is  clear  that  fault  width  controls  static  amplitude  and 
rise  time.  Beyond  a  distance  of  about  one  fault  width  from  the 
hypocenter,  the  slip  function  remains  essentially  constant  In  shape 
with  Increasing  hypocentral  distance.  This  uniformity  of  the  slip 
function  beyond  one  fault  width  Is  also  apparent  in  Figure  5  for  the 
case  w  -  4,  l  ■  16. 

In  Figure  6,  the  relations  between  fault  width  and  the  static 
offset  are  shown  In  more  detail.  Static  offset  calculated  along  the 
fault  centerline.  Is  plotted  for  the  cases  w  *  4  and  w  a  1.5  (£  *  8 
In  both  cases).  The  horizontal  lines  show  the  static  solution  for 
an  Infinitely  long  strike  slip  fault  (Knopoff,  1958).  Except  near 
the  end  of  the  fault,  the  static  slip  for  the  finite-length  faults 
Is  essentially  constant  along  the  fault  length  and  Is  very  well 
predicted  from  Knopoff s  static  solution.  At  the  center  of  the  4  x 
8,  1.5  x  8,  and  4  x  16  faults,  the  static  slip  exceeds  the  Knopoff 
solution  by  less  than  5  percent.  The  “overshoot*  phenomenon, 
observed  for  the  square  and  circular  faults,  is  not  significant  over 
most  of  the  length  of  the  long,  narrow  faults.  This  result  concurs 
with  the  numerical  results  of  Archuleta  and  Day  (1980),  which  show 
overshoot  confined  to  within  about  one  half-width  of  the  end  of  the 
fault.  Elsewhere,  the  final  shear  stress  Is  approximately  equal  to 
the  sliding  frictional  stress. 
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Slip  (m) 


Figure  3.  Computed  slip  time-histories  along  the  center-line  of  a 
rectangular  fault  with  aspect  ratio  2:1.  Scaling  is  as 
in  Figure  2. 
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Figure  5.  Computed  slip  time-histories  along  the  center-line  of  a 
rectangular  fault  with  aspect  ratio  4:1.  Slip  Is 
scaled  to  represent  a  dynamic  stress  drop  of  100  bars 
and  a  fault  length  of  16  km. 
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an  Infinitely  long  strike-slip  fault. 


Fop  a  long,  narrow  fault  model,  then,  the  static  stress  drop 
will  approximately  equal  the  dynamic  stress  drop  at,  since  little  or 
no  overshoot  occurs,  and  “undershoot"  Is  not  permitted  by  our 
model.  That  Is,  no  physical  mechanism  has  been  Incorporated  Into 
the  model  (Equation  1)  capable  of  healing  the  fault  at  a  stress 
level  higher  than  the  prescribed  sliding  frictional  level.  In 
practice,  however,  we  have  to  be  cautious  In  equating  static  and 
dynamic  stress  drops.  Seismic  estimates  of  static  stress  drop  are 
actually  estimates  of  average  static  offset  divided  by  gross  fault 
dimension.  If  an  earthquake  leaves  unbroken  patches,  or  If  some 
regions  heal  at  stress  levels  above  the  dynamic  friction  level,  then 
the  selsmlcally  Inferred  static  stress  drop  may  be  substantially 
lower  than  the  dynamic  stress  drop  at,  as  demonstrated  by  Madariaga 
(1979).  Static  stress  drop  estimates  may  constitute.  In  general,  an 
approximate  lower  bound  on  at. 

Figure  7  shows  the  relationship  between  slip  rise  time  and 
fault  width.  Rise  time  Is  plotted  along  the  fault  centerline  for  3 
rectangular  fault  calculations.  Rise  time  In  Figure  7  was  defined 
to  be  the  time  required  for  a  point  on  the  fault  to  attain  90  per¬ 
cent  of  its  final  value  of  slip.  The  horizontal  lines  represent  a 
rise  time  equal  to  the  half-width  divided  by  the  rupture  velocity. 
For  w  *  1.5,  l  •  8,  the  rise  time  at  first  decreases  with  distance 
from  the  hypocenter,  then  approaches  a  constant  level  of  about 
w/2ur.  For  w  *  4,  l  *  16,  a  constant  level  of  w/2uR  is  again 
approached  as  hypocentral  distance  incre*?**.  For  w  *"  4,  2  *  8,  the 
rise  time  again  decreases  with  distance  from  the  hypocenter,  but  the 
effects  of  the  end  of  the  fault  intervene  to  further  reduce  the  rise 
time  before  a  constant  level  can  be  clearly  established.  These 
numerical  results  predict  that  a  long,  narrow  fault  will  have  a  rise 
time  of  roughly  w/2uR  over  most  of  Its  length,  with  larger  values 
near  the  hypocenter  and  lower  values  near  the  ends. 

Actually,  In  these  simulations  the  rupture  and  shear  velo¬ 
cities  differ  only  by  about  10  percent,  so  Figure  7  could  alter¬ 
natively  be  Interpreted  as  showing  rise  time  controlled  by  1/s 
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Figure  7.  Slip  rise  time  along  the  center-line  (x  axis)  of  rectangular  faults.  Rise  time  is  defined 
as  the  time  for  slip  to  reach  90  percent  of  Its  final  vslue.  Horizontal  lines  represent 
rise  time  equal  to  the  half-width  divided  by  the  rupture  velocity. 


rather  than  l/uR.  We  Prefer  the  latter  Interpretation,  as 
explained  in  a  later  section. 

3.4  SLIP  VELOCITIES 

We  turn  our  attention  now  to  the  high  frequencies.  Here  it  is 
appropriate  to  focus  on  the  slip  velocity  function  and, 
particularly,  on  the  peak  slip  velocity.  Figure  8  shows  the  slip 
velocity  time  histories  at  selected  points  along  the  centerlines  of 
the  4  x  16  km  fault  model.  The  slip  velocities  have  been  low-pass 
filtered  to  remove  frequencies  In  excess  of  5  Hz,  which  is  close  to 
the  highest  frequency  that  can  be  reliably  computed  in  the  finite 
difference  mesh. 

Figure  8  shows  that  the  peak  slip  velocity  Initially  Increases 
with  hypocentral  distance.  However,  the  rapid  growth  in  peak  slip 
velocity  ceases  beyond  a  hypocentral  distance  of  about  one  fault 
width.  From  then  on,  the  slip  velocity  function  is  nearly  uniform 
along  the  fault  centerline.  The  figure  also  shows  slip  velocities 
observed  at  points  distributed  across  the  fault  width,  at  a  fixed 
distance  of  6  km  along  the  length  of  the  4  x  16  km  fault.  This 
figure  illustrates  the  near  uniformity  of  peak  slip  velocity  across 
the  fault  width  as  well  as  along  the  fault  length,  for  hypocentral 
distances  greater  than  about  w.  The  slight  decrease  in  peak  slip 
velocity  very  near  the  fault  edge,  evident  in  this  figure,  may  be 
due  to  the  slip  function  rise  time  near  the  edge  becoming  comparable 
to  the  rise  time  of  the  5  Hz  low-pass  filter. 

The  uniformity  of  peak  slip  velocity  is  further  illustrated  in 
Figure  9.  The  broken  curves  represent  peak  low-passed  (5  Hz)  slip 
velocities  obtained  along  the  fault  centerline  for  the  two  cases  w  » 
1.5,  l  •  8  and  w  *  4,  l  »  16  km,  respectively.  In  both  cases,  the 
peak  slip  velocity  first  Increases  rapidly  with  hypocentral 
distance,  then  quickly  settles  to  a  uniform  level  when  the 
hypocentral  distance  exceeds  w.  The  results  in  Figure  9  will  be 
further  analyzed  in  the  next  section. 
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Figure  8.  Slip  velocity  time-histories  for  the  4  x  16  km  fault. 

The  time-histories  have  been  low-pass  filtered  with  a 
5  Hz  cutoff  (corresponding  to  the  limiting  frequency 
for  which  the  numerical  method  Is  accurate).  Peak  slip 
velocity  Is  nearly  Invariant  with  position,  for 
distance  x  greater  than  fault  width. 
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Figure  9.  Peak  slip  velocity  as  a  function  of  hypocentral 
distance  along  the  fault  center-line  (x  axis).  The 
dashed  curves  represent  the  numerical  solutions  for 
1.5  x  8  km  and  4  x  16  km  rectangular  faults, 
respectively.  The  solid  curve  represents  Kostrov's 
analytical  solution  for  an  expanding  circular  crack. 
Horizontal  lines  approximate  Kostrov's  solution 
evaluated  at  radius  r  a  w.  Slip  velocities  for  the 
numerical  solutions  were  digitally  low-pass  filtered 
with  a  5  Hz  cutoff,  and  the  analytical  solution  was 
analytically  low-passed  by  convolution  with  a  boxcar 
function  of  width  l/fc  (fc  ■  5  Hz). 
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IV.  ANALYSIS  OF  THE  SLIP  FUNCTION 


4.1  THE  VELOCITY  SINGULARITY 

We  can  Interpret  the  numerical  results  for  peak  slip  velocity 
by  means  of  the  closed-form  analytic  solution  of  Kostrov  (1964)  for 
an  expanding  circular  crack.  The  analytic  solution  for  the  velocity 
discontinuity  s  on  a  circular  crack  expanding  at  a  uniform  rupture 
velocity  is 

*  at  T  +  r/\jp 

S  *  C  —  B  - - - -  H(T)  (3) 

“  VT(T+2r/UR) 

where  T  Is  the  reduced  time  (time  minus  rupture  arrival  time),  r  is 
2  2  1/2 

the  distance  (x  +  y  )  '  from  the  center  of  the  crack  to  the 
observation  point,  H  Is  the  unit  step  function,  and  C  is  a  constant 
which  equals  0.81  for  a  Poisson's  ratio  of  0.25  and  rupture  velocity 
of  0.9  s  (Dahlen,  1974).  This  solution  is  singular  at  the  rupture 
arrival  time  (except  at  the  hypocenter),  and  approaches  C^b  for  T 
large  compared  to  the  rupture  arrival  time.  In  order  to  compare 
this  solution  to  the  numerical  solutions  for  finite-width  faults,  we 
approximate  the  effect  of  a  low-pass  filter  by  averaging  the  analy¬ 
tic  solution  (Equation  3)  over  a  "cutoff"  period  l/fc.  The  re¬ 
sulting  expression  for  peak  slip  velocity  s  on  the  circular  crack  Is 

i  *  c£-  s(2rfcAjR  +  1)1/2  .  (4) 

1/2 

which,  for  fc»yr,  is  proportional  to  r  ,  That  is,  in  the 
absence  of  edge  effects  (and  nonlinearities),  peak  slip  velocity 
would  Increase  as  the  square  root  of  distance  from  the  point  of 
rupture. 

Equation  4,  with  f  equal  to  5  Hz,  is  plotted  as  a  solid 
curve  In  Figure  9.  Comparing  this  curve  with  the  peak-velocity 
curves  from  the  finite-differences  fault  simulations,  we  see  that 
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edge  effects  do  not  act  to  modify  the  peak  velocity  within  a  hypo- 

central  distance  of  approximately  one  fault  width.  Up  to  that 

distance,  the  behavior  of  peak  slip  velocity  closely  follows  that  of 

1/2 

Kostrov's  expanding  circular  crack  solution.  Increasing  as  r  . 
The  Influence  of  the  fault  edges  Is  to  terminate  this  growth  of  peak 
slip  velocity  at  a  hypocentral  distance  of  approximately  r  ■  w. 
Thus,  the  numerical  solutions  predict  that,  over  most  of  Its  length, 
a  long,  narrow  fault  will  have  a  peak  slip  velocity,  after  low-pass 
filtering  with  cutoff  frequency  fc,  of  approximately 

f  •  151 

In  (5)  we  have  assumed  fc  >>  uR/w  and  have  Introduced  the  ap¬ 

proximation  C  «ur/b,  which  Is  accurate  within  about  10  percent 
for  all  sub-shear  rupture  velocities  (Oahlen,  1974).  Equation  (5) 
Is  shown  by  horizontal  lines  on  Figure  9.  For  this  model  of 

faulting,  then,  peak  (low-passed)  slip  velocity  Is  proportional  to 
dynamic  stress-drop  at,  the  square  root  of  rupture  velocity,  and  the 
square  root  of  fault  width. 

The  above  Interpretation  of  the  numerical  solutions  permits  us 
to  "undo"  the  filtering  effect  of  the  numerical  scheme  and  char¬ 

acterize  the  slip  velocity  singularity  at  the  leading  edge  of  a 
long,  narrow  rupture.  To  simplify  the  discussion,  we  take  x»w,  so 
that  we  can  Ignore  the  rupture  front  curvature  across  the  fault 

width,  and  so  that  the  rupture  front  represents  Mode  II  (in-plane 
shear)  crack  extension.  We  expect  a  singularity  of  the  form 

s  ~*V(\jRt  -  x)"1^2  H(uRt  -  x)  (6) 

This  singular  form  is  a  universal  property  of  sub-shear-velocity 
crack  propagation  (Freund  and  Clifton,  1974;  Freund,  1979).  In 
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Equation  (6),  V  Is  the  velocity  Intensity,  and  we  can  estimate  It 
from  Equation  (5).  Recalling  that  s  Is  approximately  s  averaged 
over  time  1/fc, and  using  Equation  (6),  we  obtain  the  expression 

rj  4T 

*7  uR  (7) 

for  the  velocity  Intensity.  Thus,  the  model  predicts  a  velocity 
Intensity  proportional  to  the  square  root  of  fault  width  and  also 
proportional  to  rupture  velocity.  This  result  contrasts  with 
two-dimensional  crack  solutions.  In  which  the  velocity  Intensity  Is 
proportional  to  the  square-root  of  fault  length. 

In  addition,  we  can  use  this  result  to  estimate  the  dynamic 
stress  Intensity  factor  X  at  the  advancing  edges  of  the  fault.  The 
shear  stress  change  near  the  edge,  t  -tq,  has  the  asymptotic  form 
(Freund  and  Clifton,  1974). 

t-Tq  —  K  (x-v^t)-1^2  H(x-uRt).  (8) 


For  a  Mode  II  crack,  K  can  be  obtained  from  V  (Freund,  1979): 
2 

17?  v  • 


2ub 


3 

UR 


where  RIs  the  Rayleigh  function 
R(c) 


(9) 
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From  (7)  and  (9)  we  obtain  the  approximation 

177  *T  * 


_  a  2  R(u-) 
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(10) 
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so  that  the  dynamic  stress  Intensity  factor  Is  proportional  to  the 
square  root  of  the  fault  width  and  Is  Independent  of  the  length  of 
rupture.  This  value  of  K  happens  to  be  very  close  to  the  dynamic 
stress  Intensity  factor  for  a  two-dimensional  Inplane  shear  crack  of 
length  w;  In  fact,  the  difference  (about  8  percent  at  uR  *  0.9b) 
Is  not  significant  In  view  of  the  approximate  nature  of  our  analysis. 

Finally,  from  Equation  (10)  we  can  determine  the  energy  flux 
Into  the  propagating  rupture  front,  the  so-called  energy  release 
rate  G.  G  Is  the  energy  absorbed  per  unit  area  by  the  advancing 
crack  edge,  and  Is  given  by  (Freund,  1972) 

G  *  if—  KV  , 
from  which  we  obtain 


That  Is,  for  a  long,  narrow  fault,  the  energy  release  rate  Is 
proportional  to  the  fault  width,  rather  than  fault  length. 

4.2  THE  STATIC  SLIP 

From  Figure  6,  It  was  observed  that  the  final  slip  for  the 
rectangular  fault  model  Is  very  close  to  Knopoff's  two-dimensional 
(antiplane)  static  solution,  except  near  the  ends  of  the  fault. 
Comparing  the  two  cases  shown  in  that  figure,  it  is  evident  that  the 
length  of  the  end  region  is  proportional  to  fault  width,  rather  than 
fault  length.  The  greatest  deviation  from  the  two-dimensional 
solution  occurs  within  a  distance  w/2  of  the  ends.  We  use  these 
observations  to  construct  an  approximate  expression  which  summarizes 
the  numerical  results. 

Figure  10  shows  In  more  detail  the  numerical  results  for 
static  slip  on  the  fault  with  aspect  ratio  2.  Four  profiles  across 
the  narrow  dimension  of  the  fault  are  shown  as  dashed  curves.  The 
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uppermost  solid  curve  Is  Knopoff's  static  solution  for  the  In¬ 
finitely  long  (In  the  x  direction)  fault  of  width  w: 


*<*)  •  -r  “[i  -(£ )2  ] 

For  points  more  than  half  a  fault  width  away  from  the  end  of  the 
fault,  this  expression  Is  a  reasonable  approximation  to  the 
numerical  solution;  for  x  ■  0  and  x  ■  0.5w,  the  agreement  Is  within 
about  10  percent  over  most  of  the  fault  width.  To  modify  this 
expression  to  account  for  end  effects,  we  are  motivated  by  the  fact 
that,  very  near  the  fault  ends,  plane-strain  conditions  should  be 
approximated.  Therefore,  guided  by  the  two-dimensional  static 
solution  of  Starr  (1928)  for  a  finite.  In-plane  shear  crack,  and  the 
observation  that  the  end  region  Is  of  length  approximately  w/2,  we 
try  the  approximation 

s«  “  ‘fa) 2]  (2  -*)1/2  51/2  •  (12> 

where  £(x)  is  1  for  points  farther  than  w/2  from  the  ends,  and 
otherwise  Is  the  normalized  distance  to  the  end  of  the  fault: 

£•  U-x{  If  |f-xj  <  5 
£(x)  -  w  1  : 

U  if  U-xf>  £ 

The  approximation  (12)  is  simply  that  the  static  slip  has  the  value 
f  w  along  the  centerline,  with  an  elliptical  cross-section  across 
the  width,  multiplied  by  a  quarter-ell  ipse  taper  near  each  end.  As 
Figure  10  Indicates,  Equation  12  represents  the  numerical  solu¬ 
tions  fairly  well  over  the  whole  length  and  width  of  the  fault, 
though  somewhat  less  well  near  the  end. 

To  summarize,  we  find  that  final  slip  for  a  long,  narrow  fault 
Is  proportional  to  fault  length  rather  than  fault  width.  The  final 
slip  is  well  approximated  by  Knopoff's  solution,  except  near  the 
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ends.  The  extent  of  this  end  region  Is  also  proportional  to  fault 
width,  not  length.  Equation  12  represents  reasonably  well  the 
overall  behavior  of  the  static  offset. 

4.3  THE  SLIP  RISE  TIME 

The  approximate  expressions  deduced  above  for  the  slip 
velocity  singularity  and  static  slip  can  be  used  to  derive  an 
approximate  expression  for  the  rise  time  TR.  Integrating  the  slip 
velocity  singularity  (Equations  6  and  7)  and  equating  the  resulting 
slip  to  the  static  slip  (Equation  12)  gives  the  following  prediction 
for  the  rise  time: 

tr*  %  M^)2]12-51  «•  a3) 

assuming  x»w.  This  derivation  of  Equation  (13)  assumes  that  the 
slip  velocity  at  a  point  follows  Equation  6  until  the  static  slip 
value  Is  reached,  then  slip  terminates  abruptly.  This  is  an  over¬ 
simplification  of  the  actual  slip  function,  as  Figure  8  shows. 
However,  the  resulting  expression  for  TR  is  In  very  good  agreement 
with  the  numerical  result,  shown  In  Figure  7,  that  the  rise  time 
along  the  fault  centerline  (y*Q)  approaches  w/2'jr  for  points  more 
than  one  fault  width  away  from  the  hypocenter.  Furthermore,  the 
numerical  result  that  TR  Increases  with  decreasing  hypocentral 
distance  for  x<w,  which  is  evident  in  Figure  7,  can  be  interpreted 
In  the  same  manner.  For  x<w,  however,  Equation  (6)  should  be 
replaced  by  Equation  (3)  to  estimate  the  rise  time. 

It  Is  perhaps  surprising  that  Equation  13  involves  the  rupture 
velocity,  rather  than  the  shear-wave  velocity,  since  the  arrival  of 
shear  waves  diffracted  from  the  long  edges  of  the  fault  might  be  ex¬ 
pected  to  control  the  rise  time  at  a  point  on  the  fault.  In  fact, 
both  Day  (1979)  and  Das  (1981)  have  adequately  explained  rise  times 
for  numerical  models  of  rectangular  faults  on  the  basis  of  dif¬ 
fracted  shear-wave  arrivals,  predicting  rise  time  proportional  to 
l/s.  Since  the  shear  and  rupture  velocities  differ  by  only  10 
percent  or  so  in  our  simulations,  we  could  not  distinguish  between  a 
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1/s  and  l/uR  dependence  for  TR  on  the  basis  of  these  numerical 
results  alone.  Instead,  the  proportionality  of  rise  time  to  the 
reciprocal  of  rupture  velocity  In  Equation  13  follows  directly  from 
the  property  that  the  velocity  Intensity  at  the  leading  edge  of  the 
fault  Is  proportional  to  uR.  This  proportionality  to  uR  Is  a 
general  property  of  dynamic  cracks  running  at  sub-Rayleigh  velocity 
(see,  for  example,  Freund,  1976).  Therefore,  It  will  take  longer 
for  a  slow-running  crack  to  reach  static  equilibrium  than  for  a 
fast-running  crack  to  do  so  (recall  that  negligible  static  overshoot 
was  found  for  the  long,  narrow  faults  and  that  static  shear  stress 
Is  required  by  the  model  to  be  less  than  or  equal  to  tf).  The 
shear-wave  diffraction  effect  may  act  to  retard  the  slip  velocity 
below  that  given  by  (6)  and  (7),  but  not  to  increase  It.  For  our 
earthquake  model,  then,  an  expression  such  as  Equation  (13)  is  a 
more  appropriate  approximation  to  the  rise  time  than  would  be  ob¬ 
tained  with  s  in  place  of  uR.  The  conclusion  might  be  different 
If  some  mechanism,  such  as  velocity-dependent  friction  (Dleterich, 
1978)  were  Incorporated  Into  the  fault  model  to  permit  healing  to 
occur  at  a  stress  level  substantially  higher  than  tf. 
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V.  SUMMARY  AND  DISCUSSION 


The  approximate  behavior  of  the  slip  function  for  a  long, 
narrow  fault,  as  deduced  from  the  foregoing  analysis  of  the  numer¬ 
ical  solutions.  Is  summarized  schematically  In  Figure  11.  Near  the 
hypocenter,  the  slip  function  Initially  resembles  Kostrov's  solu¬ 
tion,  Equation  (3),  In  which  the  velocity  singularity  grows  as  the 
square  root  of  hypocentral  distance.  For  x  (distance  along  the 
fault  length)  greater  than  the  fault  width  w,  this  growth  of  the 
velocity  singularity  ceases,  and  the  slip  function  Is  then  nearly 
Invariant  with  distance  In  the  x  direction,  except  near  the  ends. 

In  this  "steady-state 11  regime,  x>w,  both  rise  time  and  static 
slip  are  proportional  to  w,  and  the  advancing  crack  edges  have 
velocity  singularities  proportional  to  Vv7.  Slip  nearly  ceases  at  a 
given  point  after  the  crack  edge  has  advanced  a  distance  of  roughly 
one  half-width  past  the  point.  The  slip-velocity  time  function  can 
be  approximated  by 

s  «/|^  &  T  ~1/2[h(T)  -  H(T-Tr)].  (14) 

In  Equation  14,  T  Is  reduced  time  t-uR/r,  and  TR  Is  the  slip 
rise  time  as  given  by  Equation  13.  This  approximation  Is  sketched 
(for  the  case  y  *  0)  In  Figure  12  along  with  the  corresponding 
shear-stress  singularity.  Also  shown  for  comparison  in  Figure  12  Is 
a  numerical  solution  for  slip  velocity  (point  E  of  Figure  8). 

The  low-frequency  characteristics  obtained  here  for  the  slip 
function,  that  is,  rise  time  and  static  slip,  are  similar  to  those 
Inferred  from  similar  calculations  by  Archuleta  and  Day  (1980)  and 
from  a  spontaneous-rupture  numerical  model  by  Das  (1981).  The  mesh 
refinement  used  to  obtain  the  numerical  solutions  in  this  study, 
however,  has  permitted  observation  of  some  previously  unresolved 
high-frequency  characteristics  of  the  solution,  as  well.  These 
high-frequency  slip  characteristics,  specifically  the  strength  of 
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11.  Sketch  summarizing  the  approximate  behavior  of  the  slip 
function  for  long,  narrow  faults.  At  a  time  greater 
than  that  required  for  the  rupture  to  cross  the  fault 
width  ( ),  slip  is  concentrated  on  two  patches,  each 
approximately  one  half-width  long,  moving  away  from  the 
hypocentef  In  opposite  directions.  The  shape  of  slip 
function  then  remains  nearly  invariant  as  these  patches 
propagate  along  the  fault  length.  In  this  steady-state 
regime,  the  static  slip  and  rise  time  are  proportional 
to  w,  and  the  slip  velocity  has  a  square  root 
singularity  with  Intensity  proportional  to  'fiT. 
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Figure  12.  Sketch  of  a  closed-form  approximation  for  the  slip 
velocity  time-history  (Equation  14)  at  y  *  0,  x  »  w, 
compared  to  a  corresponding  numerical  solution  (point  E 
of  Figure  3).  Also  sketched  Is  the  correspondl ng 
approximation  derived  for  the  shear-stress  singularity 
(Equations  8  and  9). 
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the  velocity  and  stress  singularities  at  the  leading  edge  of  the 
fault,  are  of  particular  Importance  for  understanding  the  radiation 
of  high-frequency  seismic  energy.  For  example,  Madariaga  (1977a) 
has  shown  that  "the  high  frequencies  originate  from  the  stress  and 
slip  velocity  concentrations  In  the  vicinity  of  the  fault's  edges". 
Frequencies  In  excess  of  1  Hz  are  an  Important  component  of  the 
strong  ground  motion  recorded  In  the  imnedlate  vicinity  of 
earthquakes,  and  are  also  of  Importance  for  understanding  earthquake 
ground  motion  at  regional  distances,  where  substantial  seismic 
energy  Is  observed  In  the  1-5  Hz  range. 

We  first  consider  some  Implications  of  the  dynamic  solutions 

for  kinematic  modeling  procedures  used  to  predict  ground  motion  time 

histories.  An  assumption  commonly  made  In  kinematic  modeling  of 

earthquake  ground  oration  Is  that  the  source  can  be  represented  by  a 

uniform  dislocation,  that  Is,  a  slip  function  which  Is  uniform  In 

its  amplitude  and  Its  time  dependence  over  the  entire  fault  plane 

(apart  from  a  time  delay  associated  with  the  rupture  arrival  time). 

Usually  the  slip  function  is  assumed  to  have  a  simple  ramp  time- 

history,  l.e.,  constant  slip  velocity  over  some  specified  rise 

time.  We  focus  for  now  on  the  high-frequencies  radiated  from  the 

leading  edge  of  the  fault  after  it  has  propagated  more  than  a  fault 

width  from  the  hypocenter.  In  this  case.  Equation  (14)  is  an 

appropriate  representation  of  the  slip  function  for  our  dynamic 

solutions.  The  spectral  amplitude  of  a  ramp  function  behaves 

asymptotically  at  high  frequency  as  f'1,  whereas  Equation  14 
-1/2 

implies  an  f  slip-function  spectrum  (Lighthill,  1958,  p.  52). 
Since  the  predicted  radiation  depends  linearly  on  the  assumed  slip 
function,  we  would  expect  the  ramp-function  source  representatl on  to 
be  relatively  deficient  In  Its  prediction  of  high-frequency  ground 
motion.  For  example,  one  might  attempt  to  approximate  the  slip 
function  using  a  ramp  function  In  which  the  ratio  of  static  slip  to 
rise  time  is  the  same  as  for  the  dynamic  solution  (Equation  14). 
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Then,  at  a  period  of  TR/1Q,  the  ramp  slip  function  would  be 
deficient  In  spectral  content  by  about  a  factor  of  7  relative  to 
Equation  14. 

This  spectral  comparison  between  the  ramp-function  represen¬ 
tation  and  the  dynamic  solution  may  provide  a  physical  basis  for  a 
result  obtained  by  Oel  Mar  Technical  Associates  (1979).  They 
modeled  the  1966  Parkfleld  earthquake  using  a  uniform-dislocation 
earthquake  model  with  a  ramp  slip  function,  attempting  to  fit 
spectral  characteristics  of  the  ground  motion  recorded  at  the  5 
accel erograph  stations  of  the  Chalome-Shandon  array.  They  reported 
difficulty  In  matching  observed  response  spectra  over  a  broad  period 
range  with  this  slip  function;  In  order  to  fit  recorded  short-period 
spectral  levels  of  ground  motion.  It  was  necessary  to  tolerate  a 
large  overestimate  at  long  periods.  This  result  Is  In  accord  with 
our  prediction  that  the  ramp  function  earthquake  representatl on  Is 
relatively  deficient  In  high  frequencies. 

Of  course,  a  ramp  function  suitably  scaled,  might  be  adequate 
for  modeling  ground  motion  over  a  narrow  frequency  band.  As  an 
example,  Bouchon  (1979)  successfully  synthesized  the  velocity  and 
displacement  pulses  recorded  for  the  Parkfleld  event  at  Station  2, 
using  a  uniform  slip  function  with  a  ramp  time  history.  In  this 
case,  however,  the  predominant  period  of  the  waveforms  being  modeled 
was  several  times  greater  than  the  assumed  rise  time  of  the  ramp. 

The  dynamic  solutions  reveal  a  second  difficulty  with  uniform- 
dislocation  kinematic  models,  this  one  involving  the  starting  phase 
radiated  from  the  hypocenter  when  rupture  initiates.  We  have  seen 
that  the  ramp  function  is  a  poor  representation  of  the  slip  function 
for  the  dynamic  solution  at  points  well  removed  from  the  hypocen¬ 
ter.  It  Is  tempting  to  try  to  retain  the  uniform  dislocation 
approximation,  but  alter  the  time-function  to  resemble  the  singular 
behavior  of  the  dynamic  solution  in  the  "steady-state"  regime, 
x»w.  Unfortunately,  if  the  time-function  of  a  uniform  dislocation 
model  is  chosen  to  match  the  the  dynamic  solution  at  points  far  from 
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the  hypocenter  (l.e.,  replace  the  ramp  function  by  an  expression 
such  as  (14)),  then  the  starting  phase  predicted  for  the  uniform 
dislocation  will  be  much  larger  than  that  for  the  dynamic  model.  To 
see  this,  we  use  the  expression  derived  by  Richards  (1973)  for  the 
first-motion  approximation  to  the  shear-wave  acceleration  due 
to  an  expanding  circular  crack: 

CUg 

u.  *  8  “7 - - - T?  H(t— R/®)  (15) 

'  (w|/t*  sf„2.)*  “ 

where  R  contains  the  double  couple  radiation  pattern,  e  Is  the  angle 
formed  at  the  hypocenter  by  the  fault  normal  and  the  receiver  direc¬ 
tion,  and  R  Is  the  hypocentral  distance  to  the  receiver.  We  will 
also  use  the  corresponding  expression  for  an  expanding,  uniform, 
circular  dislocation: 
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where  s(t)  is  the  assumed  slip-velocity  time  function  for  the  uni¬ 
form  dislocation.  Equation  (16)  can  be  deduced  from  Savage  (1966). 
If,  In  a  uniform-dislocation  model,  we  choose  s  in  accordance  with 

Equation  14,  Equation  16  Implies  that  the  acceleration  first-motion 
-1/2 

will  have  a  T  singularity,  whereas  the  dynamic  solution  for 
the  Initiation  phase.  Equation  15,  gives  only  a  step  discontinuity 
for  To  quantify  further  the  effect  of  the  uniform-dislocation 
approximation  In  this  case,  assuming  a  band-limited  observation  with 
cutoff  frequency  f  ,  we  substitute  Equation  14  into  Equation  16, 
average  over  a  period  l/f_,  and  take  the  ratio  of  (16)  to  (15). 
The  result  Is  that  the  uniform  dislocation  yields  a  starting  accel¬ 
eration  jghase_<wh1ch__exceeds  the  dynamic  solution  by  a  factor  of 
■/iTS,2  s.n2,)r„  fc  (we  have  Introduced  the  approximation  TR « 
w/2ur  and  uR/a  *  C).  Thus,  for  periods  much  less  than  the  rise 
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tine,  the  uniform  dislocation  model  substantially  over-predicts  the 
starting  phase  amplitudes,  relative  to  the  dynamic  solution  (except 
near  e  »  90*). 

In  summary,  the  uniform  dislocation  kinematic  models,  since 
they  do  not  Incorporate  the  spatial  variation  of  slip  present  In  the 
dynamic  solution  for  x<w,  will  either  over-predict  the  starting 
phase  or  under-predlct  the  dominant  high-frequencies  associated  with 
the  velocity  singularity.  At  least  for  the  simple  dynamic  model 
studied  here,  we  can't  construct  a  uniform  dislocation  approximation 
which  will  replicate  both  the  starting  phase  radiation  and  the 
high-frequency  behavior  of  slip  at  the  leading  edge  of  the  fault. 

The  numerical  solutions  obtained  here  also  have  tectono- 
physical  Implications.  Husselnl,  et  al_  (1975),  Oas  and  Akl  (1977b), 
and  Akl  (1979)  have  analyzed  the  role  of  fault-plane  "barriers",  or 
high-strength  segments.  In  resisting  or  arresting  rupture  growth. 
Whether  fault  growth  stops  at  a  barrier  or  rupture  continues  through 
the  barrier  depends  upon  the  barrier  strength,  as  well  as  upon  the 
stress  Intensity  factor  developed  at  the  leading  edge  of  the  fault. 
In  two-dimensional  fracture  mechanics,  the  stress  Intensity 
Increases  with  rupture  growth  as  the  square  root  of  fault  length, 
for  a  given  dynamic  stress  drop.  Therefore,  one  would  predict  that 
as  fault  length  increased,  the  likelihood  of  rupturing  barriers 
would  increase  as  well.  In  contrast,  the  three-dimensional 
solutions  demonstrate  that  the  stress  intensity  ceases  to  grow  with 
fault  length  beyond  a  distance  of  about  one  fault  width.  Equation 
(10),  derived  on  the  basis  of  the  three-dimensional  numerical 
solutions,  gives  a  stress  Intensity  factor  proportional  to  the 
square-root  of  fault  width.  The  capability  of  a  long  narrow  fault 
to  rupture  barriers  should  be  proportional  to  dynamic  stress  drop, 
and  the  square-root  of  fault  width,  but  Independent  of  rupture 
length. 

Similarly,  two-dimensional  analytical  and  numerical  solutions 
Imply  that  the  energy  release  rate  at  the  fault  edge  increases  with 
fault  length.  As  Andrews  (1976)  points  out,  those  results  would 


SYSTEMS.  SCIENCE  UNO  SOFTWARE 


/ 


require  that  longer  ruptures  be  accompanied  by  thicker  zones  of 
microcracking  or  other  Inelastic  response.  Our  three-dimensional 
results,  however.  Imply  an  energy  release  rate  proportional  to 
width,  rather  than  length  (Equation  11).  So  we  conclude  that  the 
level  of  Inelastic  response  associated  with  long,  narrow  ruptures 
will  depend  on  fault  width,  not  length. 

The  numerical  results  also  have  important  potential 
applications  to  the  problem  of  seismic  monitoring  of  nuclear 
test-limitation  treaties.  For  example,  some  of  the  methods  proposed 
for  discriminating  earthquakes  from  underground  explosions  make  use 
of  characteristic  spectral  differences  in  the  radiation  from  the  two 
types  of  events  (e.g.,  Bache,  et  a!.,  1979).  Results  such  as  these 
can  help  provide  a  theoretical  basis  for  such  discriminants. 

A  second  example  relates  to  the  seismic  estimation  of 
explosion  yield.  This  problem  requires  a  good  understanding  of 
propagation  effects  on  the  seismic  signal,  particularly  the  average 
path  attenuation  (l.e.,  t*).  Anelastlc  attenuation  is  difficult  to 
separate  observational ly  from  source  effects,  however.  In  order  to 
make  use  of  earthquake  signals  to  infer  propagation-path 
attenuation,  it  is  therefore  Important  to  have  a  good  earthquake 
source  model,  a  point  also  underscored  by  Hanks  (1981).  To  address 
these  questions  using  the  numerical  simulations,  it  will  be 
necessary  to  compute  the  radiated  seismic  signal  from  the  earthquake 
models.  The  slip  histories  obtained  here  are  sufficient  for 
synthesizing  the  radiated  waveforms,  and  this  will  be  undertaken  in 
a  subsequent  report. 
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